home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / hips / sources / 3d-tools / 3dmask.c next >
Encoding:
C/C++ Source or Header  |  1992-01-21  |  11.6 KB  |  401 lines
  1. /*
  2. % 3dmask.c --    filter an image by applying one or more masks and then
  3. %        applying another function to the various mask outputs.
  4. %        The input is in Byte, Short, Integer and Floating Point format,
  5. %        the output is floating point.
  6. %
  7. %    Copyright (c)    1990, 1991    Jin, Guojun
  8. %
  9. % usage:
  10. %    3dmask [-m filter_number] [+r] [-M] [<] input [> +o output]
  11. %            or
  12. %    3dmask [-f filter_descriptor_file] [+r] [<] input [> +o output]
  13. %
  14. %    +o output is for binary & ascii difference on PC.
  15. %    The default filter is 1.  The definition for each of these filters
  16. %    is to be found in /???/3dmasks/mask#?.@, where #? is the combinations
  17. %    of the filter_size+number_of_sets+function_number and filter type;
  18. %    @ is extension used to describe function, rotation direction, angle,
  19. %    or complicated filter name etc..
  20. %    The -f switch allows a new filter to be supplied by the user.
  21. %    The format of the filter definition file is as follows:
  22. %
  23. %    "filter name and description"
  24. %    masksize number_of-Set_of_Masks Function_name
  25. %    Number_of_masks_in_set1
  26. %    frame_position    mask--1
  27. %      .
  28. %      .
  29. %      .
  30. %    frame_position    mask--(masksize)
  31. %
  32. %    number_of_masks-(number-of-set)
  33. %    frame_position    mask--1
  34. %      .
  35. %      .
  36. %      .
  37. %    frame_position    mask--(masksize)
  38. %
  39. % where the masksize is the length of a side of all masks (which must be
  40. % square), masks are given as a sequence of integers in column-fastest order;
  41. % frame-position is for 3rd direction (frames), its range is from
  42. % -1/2 masksize to 1/2 masksize; the number-of-masks is different in each set,
  43. % since the elements of some filter frames are total zero, for fast processing,
  44. % we omit those frames, also the frame-position will missing in that(those)
  45. % set(s); the externsion is used for exactly same filter model with different
  46. % mask operations, such as Sobel, Prewitt etc.;
  47. % and the function applied to the output of the masks is chosen from:
  48. %
  49. %    1    MAXABS    - the maximum absolute value of all mask outputs
  50. %    2    MEANSQ  - the square root of the sum of the squares of all masks
  51. %    3    SUMABS  - the sum of the absolute value of all mask outputs
  52. %    4    MAX    - the maximum mask output
  53. %    5    MAXFLR    - the maximum mask output, floored at zero.
  54. %    6    MAXSFLR    - the larger of |mask-1| and |mask-2|, minus |mask-3|,
  55. %              floored at zero.
  56. %    7    MUL    - the product of the mask outputs, each floored at zero.
  57. %    8    NORM    - the first mask output normalized by the sum of the
  58. %              mask entries.
  59. %    9    DIFF    - the value of the pixel minus the normalized mask
  60. %              output. Byte format only.
  61. %    10    ORIENT    - compute orientation: 360*atan(mask2/mask1)/2*PI
  62. %
  63. % note:    important MESSAGE -- Be sure to use -O (Optimize) option to compile.
  64. %    Otherwise, you get 76% slower !!
  65. % cc -O -DMY_LIB=$(LIB_PATH) -o DESTDIR/3dmask 3dmask.c -lccs -lhipsh -lhips -lm
  66. %
  67. % AUTHOR:    Jin Guojun - LBL    12/6/90
  68. % 2/1/91:    adding short, int & float model
  69. */
  70.  
  71. #include "header.def"
  72. #include "imagedef.h"
  73. #include <math.h>
  74.  
  75. U_IMAGE    uimg;
  76.  
  77. #define    frm        uimg.frames
  78. #define    row        uimg.height
  79. #define    col        uimg.width
  80.  
  81. #ifndef    MY_LIB
  82. #define    MY_LIB        "/home/oliver/data2/users/jin/3dmasks/mask"
  83. #endif
  84.  
  85. #define    tfrm(f)        (f<0 ? 0 : (f>=frm ? frm-1 : f))
  86. #define    tcol(c)        (c<0 ? 0 : (c>=col ? col-1 : c))
  87. #define    trow(r)        (r<0 ? 0 : (r>=row ? row-1 : r))
  88.  
  89. char    *s, maskfile[128] = MY_LIB, usage[]="options\n\
  90. -m mask#    a standard mask number\n\
  91. -f maskfile    a non-standard mask file\n\
  92. +r        reverse values in all masks\n\
  93. +o output_file    for non-binary file system\n\
  94. [<] in_file [> out_file]\n";
  95.  
  96. #define    name    maskfile
  97.  
  98.  
  99. main(argc, argv)
  100. int    argc;
  101. char**    argv;
  102. {
  103. MType    fsize, i, j, k;
  104. int    mfunc,
  105.     sm, smask,    /*    # of mask sets        */
  106.     nm, nmask,    /*    # of masks in each set    */
  107.     masksz,
  108.     **mlist,    /*    mask buffers        */
  109.     mskrev=1,    /*    mask reversor    */
  110.     bp_inc, normd;    /*    normalization divisor    */
  111. int    flag,    /* out of bound flag for any one side boundary    */
  112.     boundl, boundr,
  113.     boundt, boundb,
  114.     boundf, bound,
  115.     minusd, plusd,
  116.     f=0, r, c, df, dr, dc,
  117.     framen=0,
  118.     *mval;        /*    masks working place    */
  119. float    val, *obuf, *obp;
  120. void    *buf[MaxMaskSZ];    /* maximum frames can be loaded into memory */
  121. bool    Msg=0;
  122. double    d2r = 180 / M_PI;
  123. FILE    *mfp;
  124. register int    *mskp;
  125.  
  126. format_init(&uimg, IMAGE_INIT_TYPE, HIPS, -1, *argv, "S20-1");
  127.  
  128. for (i=1, c=0; i < argc; i++, c=0)
  129.     switch(argv[i][c++]){
  130.     case '-':
  131.     switch (argv[i][c++]){
  132.     case 'M':    Msg++;    break;
  133.     case 'm':
  134.         if (!argv[i][c]){    i++; c=0;    }
  135.         strcat(maskfile, argv[i]+c);
  136.         goto    fst;
  137.     case 'f':
  138.         if (!argv[i][c]){    i++; c=0;    }
  139.         strcpy(maskfile, argv[i]+c);
  140. fst:        f = 1;    /*    filter has set up    */
  141.         break;
  142.     default:goto    info;
  143.     }
  144.     break;
  145.     case '+':
  146.     switch (argv[i][c++])
  147.     {
  148.     case 'r':    mskrev = -1;    break;
  149.     case 'o':
  150.         if (!argv[i][c]){    i++; c=0;    }
  151.         if (out_fp=freopen(argv[i]+c, "wb", stdout))    break;
  152.         message("output file %s opend error", argv[i]);
  153.     default:
  154.         goto    info;
  155.     }
  156.     break;
  157.     default:if (freopen(argv[i], "rb", in_fp) == NULL)
  158. info:        usage_n_options(usage, i, argv[i]);
  159.     }
  160.  
  161. io_test(fileno(in_fp), goto    info);
  162.  
  163. if (!f)    strcat(maskfile,"1");
  164.  
  165. (*uimg.header_handle)(HEADER_READ, &uimg, 0, 0);
  166.  
  167. if (uimg.in_form > IFMT_FLOAT)    syserr("format level should be lower than FP");
  168. uimg.o_form = IFMT_FLOAT;
  169. uimg.pxl_out = sizeof(float);
  170.  
  171. (*uimg.header_handle)(HEADER_WRITE, &uimg, argc, argv, True);
  172.  
  173. if ((mfp=fopen(maskfile,"r")) == NULL)
  174.     syserr("can't open mask file %s", maskfile);
  175. s = name;
  176. while((*s++ = getc(mfp)) != '\n');    *s = NULL;
  177. fscanf(mfp, "%d %d %d", &masksz, &smask, &mfunc);
  178.  
  179. if (masksz<1 || masksz>MaxMaskSZ || smask<1 || smask>MAXMASKSet ||
  180.     mfunc<1 || mfunc>MAXMASKFUNCS)
  181.     message("Bad filter descriptor value:    mask_size=%d, ",masksz),
  182.     syserr("sets of masks=%d, function_number=%d\n", smask, mfunc);
  183.  
  184. fsize = row * col;
  185. mlist    = (int**)zalloc((MType)smask, (MType)sizeof(*mlist));
  186. mval    = (int*)zalloc((MType)smask, (MType)sizeof(*mval));
  187. obuf    = (float*)nzalloc(fsize, (MType)sizeof (*obuf));
  188. for (i=0; i < masksz; i++)    /* locate working frame-bufffers    */
  189.     buf[i] = zalloc(fsize, (MType)sizeof(buf[0]));
  190.  
  191. for (sm=normd=0; sm < smask; sm++)    {    /* Reading Mask Sets    */
  192.     fscanf(mfp, "%d", &nm);/* get the number of masks in this set    */
  193.     mskp = (int*)nzalloc((MType)(masksz*masksz+1)*nm+1, (MType)sizeof(*mskp));
  194.     mlist[sm] = mskp;
  195.     *mskp++ = nm;        /* first value is mask numbers    */
  196.     if (nm>masksz)    syserr("%dmasks in set%d is > masksize %d",nm,sm,masksz);
  197.  
  198.     for (; nm > 0; nm--)    /* fetch mask elements    */
  199.         for (k=masksz*masksz; k>=0; k--){    /* incl frame offset    */
  200.         if (fscanf(mfp, "%d", &j) == EOF)
  201.            syserr("%s: unexpected end of the filter descriptor file\n",
  202.             argv[0]);
  203.         normd += j;    /* sum all mask elements in every sets    */
  204.         *mskp++ = j * mskrev;
  205.         }
  206. }
  207.  
  208. message("%s[%s]: using filter:    %s\n", *argv, Mversion, name);
  209.  
  210. plusd    = masksz / 2;        /* set central position    */
  211. minusd    = plusd - masksz + 1;
  212. boundl=boundt=boundf = -minusd;    /* compute boundaries of where the    */
  213. boundr    = col - plusd - 1;    /* mask overlaps the image completely    */
  214. boundb    = row - plusd - 1;    /* for more efficient convolution &    */
  215. bound    = frm - plusd - 1;    /* calculate bufp move difference for    */
  216. bp_inc    = col - masksz;        /* down one pos where not at boundary.    */
  217.  
  218. /*************************************************
  219. ********    start    masking        *********/
  220.  
  221. for (df=0; df<plusd; df++)    /* get first half set of frames    */
  222.     if ((f=upread(buf[df], uimg.pxl_in, fsize, in_fp)) != fsize)
  223.     syserr("unexpected end of image_file");
  224. for (f=0; f<frm; f++)
  225. {
  226.     obp = obuf;
  227.     if (f+df < frm)    /* get next frame    */
  228.        if (upread(buf[(f+df) % masksz], uimg.pxl_in, fsize, in_fp) != fsize)
  229.         syserr("unexpected end_of_file frame{%d}", f);
  230. #ifdef    _DEBUG_
  231.     else if(Msg)    message("masking frame[%d], ef[%d], lc[%d]\n",
  232.             f, f+df, (f+df) % masksz);
  233.     message("now: buf=%d, tfrm=%d\n",
  234.         buf[tfrm(f+framen)%masksz], tfrm(f+framen));
  235. #endif
  236.     for (r=0;r < row; r++)
  237.         for (c=0; c < col; c++){
  238.         flag = (r < boundt || r > boundb ||
  239.             c < boundl || c > boundr ||
  240.             f < boundf || f > bound);
  241.  
  242.         for (sm=0; sm < smask; sm++)
  243.         {
  244.             k = 0;
  245.             mskp = mlist[sm];
  246.             nmask = *mskp++;
  247.             for (nm=0; nm < nmask; nm++){
  248.             framen = *mskp++;
  249.             switch (uimg.in_form){
  250.             case IFMT_BYTE:{
  251.             register byte*    bufp;
  252.                 if (flag) {
  253.                 byte*    tmpp = buf[tfrm(f+framen) % masksz];
  254.                 for (dr = minusd; dr <= plusd; dr++)
  255.                 {
  256.                 bufp = tmpp + trow(r+dr) * col;
  257.                 for (dc=minusd; dc <= plusd; dc++)
  258.                   k += *mskp++ * (*(bufp + tcol(c+dc)) & 0xFF);
  259.                 }
  260.                 }
  261.                 else{
  262.                 bufp = (byte*)buf[(f+framen) % masksz] + (r + minusd)*col + c + minusd;
  263.                 for (dr=minusd; dr <= plusd; dr++)
  264.                     {    for (dc=minusd; dc <= plusd; dc++)
  265.                         k += *mskp++ * (*bufp++ & 0xFF);
  266.                     bufp += bp_inc;
  267.                     }
  268.                 }
  269.             }
  270.             break;
  271.             case IFMT_SHORT:{
  272.             register short*    bufp;
  273.                 if (flag) {
  274.                 short*    tmpp = buf[tfrm(f+framen) % masksz];
  275.                 for (dr = minusd; dr <= plusd; dr++)
  276.                 {
  277.                 bufp = tmpp + trow(r+dr) * col;
  278.                 for (dc=minusd; dc <= plusd; dc++)
  279.                   k += *mskp++ * *(bufp + tcol(c+dc));
  280.                 }
  281.                 }
  282.                 else{
  283.                 bufp = (short*)buf[(f+framen) % masksz] + (r + minusd)*col + c + minusd;
  284.                 for (dr=minusd; dr <= plusd; dr++)
  285.                     {    for (dc=minusd; dc <= plusd; dc++)
  286.                         k += *mskp++ * *bufp++;
  287.                     bufp += bp_inc;
  288.                     }
  289.                 }
  290.             }
  291.             break;
  292.             case IFMT_LONG:{
  293.             register int*    bufp;
  294.                 if (flag) {
  295.                 int*    tmpp = buf[tfrm(f+framen) % masksz];
  296.                 for (dr = minusd; dr <= plusd; dr++)
  297.                 {
  298.                 bufp = tmpp + trow(r+dr) * col;
  299.                 for (dc=minusd; dc <= plusd; dc++)
  300.                   k += *mskp++ * *(bufp + tcol(c+dc));
  301.                 }
  302.                 }
  303.                 else{
  304.                 bufp = (int*)buf[(f+framen) % masksz] + (r + minusd)*col + c + minusd;
  305.                 for (dr=minusd; dr <= plusd; dr++)
  306.                     {    for (dc=minusd; dc <= plusd; dc++)
  307.                         k += *mskp++ * *bufp++;
  308.                     bufp += bp_inc;
  309.                     }
  310.                 }
  311.             }
  312.             break;
  313.             case IFMT_FLOAT:{
  314.             register float*    bufp;
  315.                 if (flag) {
  316.                 float*    tmpp = buf[tfrm(f+framen) % masksz];
  317.                 for (dr = minusd; dr <= plusd; dr++)
  318.                 {
  319.                 bufp = tmpp + trow(r+dr) * col;
  320.                 for (dc=minusd; dc <= plusd; dc++)
  321.                   k += *mskp++ * *(bufp + tcol(c+dc));
  322.                 }
  323.                 }
  324.                 else{
  325.                 bufp = (float*)buf[(f+framen) % masksz] + (r + minusd)*col + c + minusd;
  326.                 for (dr=minusd; dr <= plusd; dr++)
  327.                     {    for (dc=minusd; dc <= plusd; dc++)
  328.                         k += *mskp++ * *bufp++;
  329.                     bufp += bp_inc;
  330.                     }
  331.                 }
  332.             }
  333.             }    /* end switch (uimg.in_form)    */
  334.             }    /* end for nm    */
  335.             mval[sm] = k; /* / (masksz * masksz);*/
  336.         /* need not divid by size since must use scale_gray anyway */
  337.         }    /* end for mask    set    */
  338.  
  339.         switch(mfunc) {
  340.         case MASKFUNC_MAXABS:
  341.             k = abs(mval[0]);
  342.             for (i=1; i < smask; i++)
  343.                 if(abs(mval[i]) > k)    k = abs(mval[i]);
  344.             val = k;
  345.             break;
  346.         case MASKFUNC_MEANSQ:
  347.             for (i=k=0; i < smask; i++)
  348.                 k += mval[i]*mval[i];
  349.             val = sqrt((double)k);
  350.             break;
  351.         case MASKFUNC_SUMABS:
  352.             for (i=k=0; i < smask; i++)
  353.                 k += abs(mval[i]);
  354.             val = k;
  355.             break;
  356.         case MASKFUNC_MAXX:
  357.         case MASKFUNC_MAXFLR:
  358.             k = mval[0];
  359.             for (i=1; i < smask; i++)
  360.                 if (mval[i] > k)    k = mval[i];
  361.             val = (k>0 || mfunc == MASKFUNC_MAXX) ? k : 0;
  362.             break;
  363.         case MASKFUNC_MAXSFLR:
  364.             k = abs(mval[0]) > abs(mval[1])
  365.                 ? abs(mval[0]) : abs(mval[1]);
  366.             k -= abs(mval[2]);
  367.             val = k>0 ? k : 0;
  368.             break;
  369.         case MASKFUNC_MUL:
  370.             k = 1;
  371.             for (i=0; i < smask; i++)
  372.                 k *= mval[i]>0 ? mval[i] : 0;
  373.             val = k;
  374.             break;
  375.         case MASKFUNC_NORM:
  376.             val = ((float)mval[0])/normd;
  377.             break;
  378.         case MASKFUNC_DIFF:
  379.             val = *((byte*)buf[f % masksz]+uimg.pxl_in*(r*col+c)) -
  380.                 (float)mval[0] / normd;
  381.             break;
  382.         case MASKFUNC_ORIENT:
  383.             if (mval[0] || mval[1])
  384.                 val = d2r * atan2( (double)mval[1],
  385.                     (double)mval[2] );
  386.             else    val = 0;
  387.             break;
  388.         default:
  389.             syserr("bad function type?%d", mfunc);
  390.         }
  391.         *obp++ = val;
  392.         }    /* end for coln    & end for row    */
  393.  
  394.     if (fwrite(obuf, sizeof(*obuf), fsize, out_fp) != fsize)
  395.         syserr("<#%d>: write error", f);
  396. #ifdef    _DEBUG_
  397.     else if(Msg)    message("frame %d    done\n", f);
  398. #endif
  399. }    /*    end for frame    */
  400. }
  401.